Collisionless relaxation in gravitational systems: Prom violent relaxation to 

gravothermal collapse 
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Theory and simulations are used to study collisionless relaxation of a gravitational iV-body system. 
It is shown that when the initial one particle distribution function satisfies the virial condition - 
potential energy is minus twice the kinetic energy - the system quickly relaxes to a metastable state 
described quantitatively by the Lynden-Bell distribution with a cutoff. If the initial distribution 
function does not meet the virial requirement, the system undergoes violent oscillations, resulting in 
a partial evaporation of mass. The leftover particles phase separate into a core-halo structure. The 
theory presented allows us to quantitatively predict the amount and the distribution of mass left in 
the central core, without any adjustable parameters. On a longer time scale tg ~ N collisionless 
relaxation leads to a gravothermal collapse. 



Since the pioneering works of Boltzmann and Gibbs, 
systems with long range interactions have been a ma- 
jor stumbling block to the development of statistical me- 
chanics 1]. The difficulty was already well appreciated 
by Gibbs, who has noted that the equivalence between 
statistical ensembles breaks down when the interparti- 
cle potentials decay with exponents smaller than the 
dimensionality of the embedding space Q. When this 
happens, systems exhibit some very unusual properties 
which appear to violate the second law of thermodynam- 
ics. For example, confined non-neutral plasmas are found 
to phase separate into coexisting phases of different tem- 
peratures 0], while the self-gravitating systems, such as 
elliptical galaxies, are characterized by a negative specific 
heat [|J]. The explanation for these counterintuitive re- 
sults lies in the fact that when the interactions are long 
ranged, thermodynamic equilibrium is never reached [5j 
and the laws of equilibrium thermodynamics do not ap- 
ply- 

In the limit in which the number of particle goes to 
infinity (N — > oo), while the total mass and charge are 
kept fixed — the so called thermodynamic limit — the 
collision duration time diverges, and the dynamical evo- 
lution of non-neutral plasmas and gravitational systems 
is governed exactly by the collisionless Boltzmann — or 
as it is known in plasma physics, Vlasov equation Q. 
This equation never reaches a stationary state — the spa- 
tiotemporal evolution continues ad infinitum on smaller 
and smaller length scales, while the one particle posi- 
tion and velocity distribution function evolves in time 
as an incompressible fluid. In practice, however, since 
there is always a minimum resolution most systems do 
appear to evolve to a well defined stationary state. This 
state, however, is very different from the normal ther- 
modynamic equilibrium characterized by the Maxwell- 
Boltzmann distribution — it explicitly depends on the 
initial distribution of the particle positions and veloci- 
ties. 

Forty years ago Q, Lynden-Bell argued that al- 
though the fine-grained distribution function of positions 
r and velocities v, /(i, r, v), never reaches equilibrium, 
the coarse-grained distribution function /(/;, r, v, ), av- 



eraged on microscopic length scales, relaxes to a meta- 
cquilibrium with /(r, v). Since in practice the very small 
length scales can not be resolved experimentally, obser- 
vations and simulations can only provide us with the 
information about the coarse-grained distribution func- 
tion /(r,v). To obtain /(r, v) we divide the phase 
space into macrocells of volume d 3 rd 3 v, which are in 
turn subdivided into v microcells, each of volume h 3 . 
The initial distribution function /o(r, v) is discretized 
into a set of levels rjj, with j = I... p. The incom- 
pressibility of the Vlasov dynamics requires that at any 
time t each microcell contains at most one discretized 
level rjj and that the overall hypervolume of each level 
liVj) = I <5(/(*) r ) v ) — ?7j)d 3 rd 3 v, be preserved by the 
dynamics. We denote the fraction of the volume of the 
macrocell at (r, v) occupied by the level j aspw(r, v). Us- 
ing a standard combinatorial procedure [H, 0, [|| it is now 
possible to associate a coarse-grained entropy S with the 
distribution of {pj}. Lynden-Bell argued that the colli- 
sionless relaxation should lead to the density distribution 
of levels which is most likely, i.e. the one that maximizes 
the coarse-grained entropy, consistent with the conserva- 
tion of all the dynamical invariants — energy, momen- 
tum, angular momentum, and the hypervolumes ^{rjj). 
In terms of the volume fractions {pj}, the stationary dis- 
tribution function becomes /(r, v) = Y^,j r)jpj{r, v). If 
the initial distribution has only one level p = 1 (is water- 
bag), 

/o(r,v) = ?7i0(r m - r)Q(v m - v) (1) 
where Q(x) is the Heaviside step function and rji = 
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the maximization procedure is partic- 



ularly simple, yielding a Fermi-Dirac distribution, 

/(r, v) = mP (r, v) = e/3[6(r|V) _ M] - 1 • (2) 

In the expression above, e is the mean energy of particles 
at position r with velocity v. (3 and p, are two Lagrange 
multipliers required by the conservations of energy and 
the number of particles, 



d 3 rd 3 v e(r,v)/(r,v) = e , 



(3) 
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J d 3 rd 3 v/(r,v) = 1 , 

where eo is the energy per particle of the initial distribu- 
tion and the units are such that h = 1. By analogy with 
the usual Fermi-Dirac statistics, we define the effective 
temperature of a stationary state T as (3 — 1/fcgT. This 
temperature should not be confused with the standard 
definition of temperature in terms of the average kinetic 
energy - the latter being valid only for classical systems 
in thermodynamic equilibrium. In the thermodynamic 
limit, the gravitational potential f> of N particles with 
the total mass M satisfies the Poisson equation 

V 2 = 4?rGmn(r), (4) 

were m — M/N and n(r) = N J f d 3 v is the parti- 
cle number density. The Poisson equation (fj| and the 
equations (|2l3p form the basis of Lynden-BelPs violent 
relaxation theory 0, H, [l(| • The idea is that the origi- 
nal distribution f — which is far from equilibrium i.e. 
is statistically unlikely — will relax rapidly to /(r, v), 
thus maximizing the coarse grained entropy. In prac- 
tice, however, what is found is that self-gravitating sys- 
tems usually relax to structures characterized by dense 
cores surrounded by dilute halos, the distribution func- 
tions of which are quite different from Lynden-BelPs /. 
The failure of the theory was attributed to the fact that 
the violent relaxation occurs on very fast dynamical time 
scale and the system does not have time to explore all 
of the phase space to find the most probable configu- 
rations pH . Recent work on non- neutral plasmas [3], 
however, provides a very different picture, ft has been 
found that confined non-neutral plasmas also relax to 
a core-halo structure. In that case, however, the halo 
production has been clearly shown to be the result of 
parametric resonances arising from the macroscopic bulk 
oscillations [l2|. If the initial distribution is constructed 
in such a way as to suppress macroscopic oscillations, the 
resulting stationary state was found to be precisely the 
one predicted by the Lynden-Bell theory [3j. It is rea- 
sonable, therefore, to suppose that a similar mechanism 
will be at work for the self-gravitating systems as well. 
Strong oscillations will lead to parametric resonances — 
a form of a non-linear Landau damping [l3j — which will 
transfer a large amount of energy to some particles at the 
expense of the rest. These particles will either escape to 
infinity (evaporate) or will form a dilute halo which will 
surround the central core. 

To test this theory, we first consider the case in which 
the macroscopic oscillations are suppressed. This can 
be achieved by forcing the original distribution to sat- 
isfy the virial condition 2K = — U, the virial number 
1Z = —2K/U is one, where K is the total kinetic en- 
ergy and U is the total gravitational energy. To simplify 
the discussion, we will restrict our attention to the initial 
distributions of the water-bag form (p = I). For these 
distributions the virial condition reduces to the require- 
ment that v m = \J GM/r m and the average energy per 



particle is eo = jQ mv m ~ f 7 Mm - We expect that under 
these conditions /o will relax to the distribution given 
by Eq. |J2J), with e(r, v) = mv 2 /2 + m<j){r), subject to 
constraints of Eqs. (J3]). There is, however, one difficulty. 
Since the gravitational potential decays to zero at large 
distances, Eq.Q requires that at any finite temperature 
there should be a non-vanishing particle density over all 
space. The normalization conditions ([3]), therefore, can- 
not be satisfied in an infinite space. However, if we con- 
fine our attention to sufficiently short times, before a sig- 
nificant number of particles has a chance to escape from 
the main clusters — in practice this time is very large 
when the virial condition is satisfied — the normaliza- 
tion problem can be avoided by artificially restricting the 
particle positions to lie within a sphere of radius R. The 
situation here is very similar to the one encountered in 
the theory of electrolyte solutions [1J|. In that case, the 
canonical partition function of an ionic cluster is found to 
diverge and a cutoff has to be introduced to obtain finite 
results. The divergence is a natural consequence of the 
fact that at any finite temperature ionic clusters are un- 
stable and will fall apart after a sufficiently long time. On 
short time scales, however, the dynamics of ionic clusters 
is well described by a statistical theory with a cutoff. Fur- 
thermore, the thermodynamics of electrolyte solutions at 
low temperatures is found to be completely insensitive 
to the precise value of the cutoff used [l5| . We find the 
same is true for the gravitational systems as well. In 
the infinite time limit, a gravitational cluster satisfying a 
virial condition is unstable and some particles will slowly 
evaporate. On "short" time scales, however, the cluster 
properties are well described by a statistical theory with 
a cutoff. The precise value of the cutoff is unimportant 
— as long as it is not too large. In our calculations we 
have taken the cutoff to be at R — 10r m , but all the re- 
sults remain visibly unaffected if we replace this by 5r™ 
or 100r m . The cutoff- Lynden-Bell distribution (cLB) 
is then given by /cL_b(i", v) = /(r J v)0(i? — r). It is also 
possible to use an energy cutoff [16| , but for the purposes 
of the present calculation this is not necessary. We now 
iteratively solve the Poisson equation Eq. ^ with the 
distribution f C LB subject to the conservation equations 
([3]). Integrating the Fermi-Dirac distribution over all ve- 
locities and taking advantage of the radial symmetry of 
the distribution JTJ) the Poisson's equation (01 takes the 
form 



(5) 

where (3 = /3m, f = fj,/m and Li n (x) is the n th polylog- 
arithm function of x. This nonlinear equation is solved 
numerically with the boundary conditions (j>'(r = 0) = 
and <p(r — ► oo) = 0. The gravitational potential <p(r) de- 
pends parametrically on /i and ft, which are determined 
using the conservation equations ([3]). Integrating over 
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the velocities these become 

^l L i 3/2 (-e^-^ r*dr = e , (6) 

16 ^n^ £ r 2 L l3/2 (-e^-^)dr = 1. 

To compare the theory with the simulations, we cal- 
culate the number particles inside shells located between 
r and r + dr, N(r)dr = AnNr 2 dr J d 3 v/(r, v). In the 
simulations 20, 000 particles were initially distributed ac- 
cording to the water-bag distribution Eq. ([I]) and then 
allowed to evolve in an infinite space in accordance with 
Newton's equations of motion. To avoid the collisional 
effects and to speed up the simulations, the forces were 
calculated using the mean-field Gauss law. As discussed 
earlier, this procedure becomes exact in the thermody- 
namic limit. In Fig. 1 the solid lines are the values 
of N(r)/N obtained using the theory above, while the 
points are the results of the dynamics simulation, the 
distances are measured in units of r m and the dynamical 
time scale is tjj = y/r^JGM. An excellent agreement 
is found between the theory and the simulations, with- 
out any adjustable parameters. To further explore the 
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the main cluster. Following the violent relaxation, \ be- 
gins to increase again. The rate of this increase depends 
strongly on the number of particles in the system, Fig 
2a. We define tg as the time at which the (violently) re- 
laxed distribution begins to deviate from the cLB form by 
1%, x = 0.01. This time depends on TV as t g sa 4Nt d . 
Scaling the simulation time with tg gives an excellent 
collapse of the data on a single universal curve, see Fig. 
2b. The fact that tq diverges with N implies that in 



0.04 




f/lOOOx 



^#=1000 


(b) - 


N= 2000 




-o- N = 4000 




-^N = 20000 









0.0 1.0 2.0 3.0 4.0 5.0 

t/x G 



FIG. 1: Mass distribution for 1Z = 1: solid curve is obtained 
using the cLB distribution and the points are the result of 
dynamics simulation. There is no halo, all mass is in the 
core. 

dynamics of the relaxation process, we define the tempo- 
ral deviation of the density distribution N(r, t) from the 
stationary cLB value, N c lb{t), 

X® = jjs J o [N(r,t)-N cLB (r)}\ (?) 

The inset of Fig. 2a shows that after a very short time in- 
terval of tr « to, the original distribution /o quickly 
relaxes to the cLB form. Furthermore the relaxation 
time tr is independent on the number of particles in 
the system — this is precisely the Lynden-Bell's violent 
relaxation regime. The metastable cLB distribution per- 
sists until a finite fraction of particles evaporates from 



FIG. 2: (a) x(t) f° r different number of particles in the system. 
Inset shows the violent relaxation regime which occurs on a 
very short time scale td and is independent of N. Following 
the violent relaxation, the system undergoes small oscillations 
which die out after about IOOtd . At this time the distribution 
is precisely of the cLB form. On a longer time scale tq the 
system undergoes a gravothermal collapse. The rate of this 
collapse is a linear function of N. (b) When the time is scaled 
with tg, all the data in (a) collapses on one universal curve. 

the thermodynamic limit the cLB distribution will last 
forever! We conclude that if the virial condition is satis- 
fied and the macroscopic oscillations are s upp ressed, the 
phase-mixing — linear Landau damping [17j ] — mecha- 
nism is extremely efficient to produce a local ergodicity. 
For times larger than tg , the phase space incompressibil- 
ity condition pi (r, v) < 1 is violated and the system un- 
dergoes a slow gravothermal collapse. We note, however, 
that since the time scale tg, is larger than the Chan- 
drasekhar time tch ~ Td-/V "/ 'ln(iV), the binary collisions 
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omitted in our simulations must be explicitly taken into 
account to study this regime. 

Small deviations from the virial condition 0.8 < 7Z < 
1.2 result only in weak oscillation which are not sufficient 
to produce significant parametric resonances. Thus, we 
find that for this range of virial numbers the cLB theory 
remains in good agreement with the dynamics simula- 
tions. For larger deviations from 1Z = 1, the situation 
changes dramatically. In these cases, the initial distribu- 
tion undergoes violent oscillations resulting in a partial 
mass evaporation and a halo production. A fraction of 
the particles quickly gain enough energy from the reso- 
nances to completely escape from the main cluster (evap- 
orate), while the other fraction gains only enough energy 
to move away from the core, remaining gravitationally 
bound to it. This latter class forms a dilute halo sur- 
rounding the dense central core. The evaporation and 
halo production progressively cool down the cluster until 
all the collective oscillation cease at T = 0. The parti- 
cles left in the core should then be in the energy ground 
state, with their distribution function given by that of a 
fully degenerate Fermi gas / c (r, v) = r/i 8(/i — e). This 
is precisely what is found when the cutoff R in the cLB 
distribution is extended to infinity. In this limit the cLB 
distribution splits into two domains — a compact zero 
temperature core described by / c (r, v) plus an evapo- 
rated fraction of zero energy particles at infinity. Inte- 
grating over velocities, the Poisson equation becomes 

1 d 2 d<j> 32GMir 2 V2vi,~ ,,3/2 . » 

^~dr T = 3 {jl ~ 0W) 6(M " ^ (r)) ' 

where the value of fi is determined by the energy conser- 
vation, 



The work in this direction is now in progress. On a time 
scale larger than tq the core is, once again, found to 
undergo a gravothermal collapse. 

Traditionally the failure of the Lynden-BelPs theory 
was attributed to the fact that for gravitational systems 
relaxation occurs on a short dynamical time scale rr> so 
that the system has no time to explore all of the phase 
space to find the most "probable" configuration. Present 
work, however, provides a very different picture. Strong 
oscillations lead to propagating density waves [HI and 
to parametric resonances which force some particles into 
statistically improbable regions of the phase space. These 
regions do not mix with the rest of the system. When 
the oscillations (parametric resonances) are suppressed, 




FIG. 3: Mass distribution for TZ = 1.9 - points are the result 
of the simulation. Solid curve obtained using f c gives the 
mass distribution inside the core. The theory predicts that 
54% of the mass will be in the halo or will evaporate, which is 
in agreement with the simulations. There are no adjustable 
parameters 
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r 2 (fL-cl>(r)f/ 2 G(il-<l)(r))dr = e . 



(9) 

The norm of f c then gives the amount of mass left in the 
central core after the process of collisionless relaxation is 
completed. Fig. 4 shows that the theoretically predicted 
f c is in excellent agreement with the core mass distribu- 
tion obtained using the dynamics simulations. However, 
in order to have a complete account of the halo mass dis- 
tribution a more detailed dynamical study is necessary. 



the mixing is very efficient and the predictions of the 
Lynden-Bell theory are verified quantitatively. Outside 
the virial condition, strong oscillations lead to a par- 
tial mass evaporation and a core-halo coexistence. The 
theory presented here gives a quantitatively accurate ac- 
count of the core, it also allows us to predict the amount 
of mass which will be lost to the halo production and 
evaporation. 
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